Imputation of the data

It is highly recommended to “impute” the genotype likelihood files prior to running the gwas. What does impute mean?

Imputation infers missing genotypes from sampled sites and has been shown to improve the accuracy of Another way to read this is that GWAS is sensitive to missing data. So a good double target approach would be to determine how much missing data we have….and remove samples that are mostly missing (>50%) and impute the rest. In this tutorial we will be removing the missing data and then imputing the remaining samples.

How do you know whether you have missing data or not?

We can use this command here which uses bcftools to first extract the sample names, and genotypes, and then uses awk to sum the genotypes and divide by that sum which returns the proportion of missing data for each sample. You will see that we have a number of samples with missingness greater than 50%!

paste <(bcftools query -f '[%SAMPLE\t]\n' 5.n190_moreSNPS_gwas_thinned.recode.vcf | head -1 | tr '\t' '\n') <(bcftools query -f '[%GT\t]\n' 5.n190_moreSNPS_gwas_thinned.recode.vcf | awk -v OFS="\t" '{for (i=1;i<=NF;i++) if ($i == "./.") sum[i]+=1 } END {for (i in sum) print i, sum[i] / NR }' | sort -k1,1n | cut -f 2)

You can find more helpful scripts like the one above here

Filter based on missingness

We can see that there are a few samples at the end there that have a pretty high proportion of missing data. Let’s remove those from our analyses and impute the rest.

Let’s remove individuals that are more than 50% missing data.

bcftools view -S ^<(paste <(bcftools query -f '[%SAMPLE\t]\n' 5.n190_moreSNPS_gwas_thinned.recode.vcf | head -1 | tr '\t' '\n') <(bcftools query -f '[%GT\t]\n' 5.n190_moreSNPS_gwas_thinned.recode.vcf| awk -v OFS="\t" '{for (i=1;i<=NF;i++) if ($i == "./.") sum[i]+=1 } END {for (i in sum) print i, sum[i] / NR }' | sort -k1,1n | cut -f 2) | awk '{ if ($2 > 0.50) print $1 }') 5.n190_moreSNPS_gwas_thinned.recode.vcf| gzip > 5.n190_moreSNPs_gwas_thinned_nomissing.vcf.gz

Impute with beagle

Now we will impute the remaining samples that have less missing data (some had very little missing data). This takes less than 4 minutes. If we’re short on time we may skip to the next step. Most likely you would run this on a larger data set and it would thus take a little longer. For example, in a test run on a dataset with >500k SNPs it took ~20 minutes to run. You could let it run on your screen like this, but I put it in a batch script.

wget https://faculty.washington.edu/browning/beagle/beagle.27Jan18.7e1.jar

java -jar beagle.27Jan18.7e1.jar gt=5.n190_moreSNPs_gwas_thinned_nomissing.vcf out=imputed_thinned_nomissing

Format data with plink

Gemma requires our data to be in a bimbam format so we will convert it to the necessary format with plink.

Convert the vcf file to .bed format with plink

plink --make-bed --vcf imputed_thinned_nomissing --chr-set 95 --allow-extra-chr

This produces four files. A .bed, .bim, .fam, and a log file

The .bim file is a binary so we can’t look at it but let’s take a look at the .fam file. This file has a column for ‘phenotype’ but there wasn’t phenotypic information in our vcf file, so that column shows -9 for missing data.

We now need to add our phenotypes to the .fam file. We will use the awk command below to achieve this.

awk 'FNR==NR{a[NR]=$1;next}{$6=a[FNR]}1' pheno_bin plink.fam > test.fam

head *.fam

You should see a column of -9 replaced by 1s or 0s (you’ll see the 0s if you use tail)

# then replace our original fam with the test.fam

mv test.fam plink.fam

install packages in R

Navigate to the Open OnDemand app to access RStudio here. You will need to login with your northeastern account information.

We’ll be using the package qqman to make a manhattan plot.

install.packages("qqman")

## Error in contrib.url(repos, "source"): trying to use CRAN without setting a mirror

library(qqman)

make a manhattan plot

Our GWAS results are contained within the file labeled: gwas_kinship.assoc.txt. We’ll start by reading that into R.

#change our directory to our scratch folder

#setwd("/scratch/s.caplins/training_gwas/")

gwas_df<-read.table("/Users/s.caplins/Documents/Bootcamp_training_materials/training_gwas/output/gwas_kinship.assoc.txt", header = T)

Take a look at the dataframe. The test results are located in the last three columns. We have results from a wald test p_wald, a likelihood ratio test p_lrt, and a score test p_score

head(gwas_df)

##   chr rs    ps n_miss allele1 allele0    af        beta         se   logl_H1
## 1   1  .  1360      0       G       T 0.191 -0.01038153 0.08048139 -31.08002
## 2  12  .  1762      0       T       C 0.086  0.03106614 0.11176640 -31.04891
## 3  13  .   266      0       A       C 0.026  0.17880600 0.20707290 -30.70759
## 4  13  .  3748      0       T       A 0.020  0.14090390 0.24282750 -30.91605
## 5  13  .  9158      0       A       C 0.013  0.19685070 0.28295870 -30.84084
## 6  13  . 11658      0       T       C 0.039  0.20810010 0.16631610 -30.29300
##   l_remle l_mle    p_wald     p_lrt   p_score
## 1   1e+05 1e+05 0.8977132 0.8959991 0.8963596
## 2   1e+05 1e+05 0.7818213 0.7782398 0.7790819
## 3   1e+05 1e+05 0.3906564 0.3827219 0.3867291
## 4   1e+05 1e+05 0.5634997 0.5569449 0.5591778
## 5   1e+05 1e+05 0.4888041 0.4815096 0.4844282
## 6   1e+05 1e+05 0.2147909 0.2071657 0.2134913

hist(gwas_df$p_lrt)

Everything looks good to proceed to making our manhattan plot.

the function manhattan requires each SNP to have it’s own “name”. Let’s make a vector for rownumbers that start with the letter r. And we should do some correction for multiple testing, though you will see we don’t find much with this little dataset.

gwas_df$SNP<-paste("r",1:length(gwas_df$chr), sep="")

gwas_df$pvalue<-pchisq(gwas_df$p_lrt, df=1, lower=F)

padj <- p.adjust(gwas_df$pvalue, method = "BH")
alpha <- 0.05
outliers <- which(padj < alpha)
length(outliers)

## [1] 0

hist(padj)

manhattan(gwas_df, chr="chr", bp="ps", p="p_lrt")

Let’s look at a qq-plot of our pvalues to check the model fit

qqnorm(gwas_df$p_score)